  function ocell = out_moments_both( G1 , HH , SDX , ZZ ,CC,  sercell , shocell, in , Y , info ); 
% ==================================================================================== 
% OUT_MOMENTS_BOTH 
% function ocell = out_moments_both( G1 , HH , SDX , ZZ , sercell , shocell, Y , in ); 
% 
% Combines the evaluation of moments from three sources 
% 1) DATA 
% 2) SIMULATION 
% 3) MODEL 
% 
% 1) Computed only if Y input (data vector) provided 
% Code combines OUT_MOMENTS.m and OUT_MOMENTS_SIM.m 
% Allows use of same functions although does not permit BP filtered data 
% Hence, will be used only for the observable ZZ matrix used in the
% estimation 
% Inputs 
% ======
% gensys solution to the system in the form 
% y(t) = ZZ*s(t)
% s(t) = G1*s(t-1) + HH*e(t)  V( e(t) ) = SDX'*SDX 
% 
% COMPUTE MOMENTS OF THE SOLUTION FROM GENSYS 
% 
% Compute the std, vcv, correlation, ac(1) and variance decomposition
% 
% If data vector is provided compute the equivalent data moments 
% except for the variance decomposition 
% 
% Output is OCELL 
% 
% Alejandro Justiniano 6/19/2005 
% 10/10/2006    Variance Decompositions at Different Horizons 
% ======================================================================
nper=10; 
cucd = cd; 
% Initial Checks 
% ---------------
nz = size(ZZ,1); 
nx = size(SDX,1); 
temp = max( abs( eig(G1) ) ); 
if temp > 0.98  
    warning('Warning! Root of G1 is above 0.98 numerical problems may occur'); 
end 
sercell = sercell(:);
if length( sercell ) ~= nz 
    error('Mistmatch length SERCELL and zmat matrix')
end 
if isempty( shocell ) == 1;
    shocell = fnumcell( 'Shock ', 1:nx ) ;
    shocell = shocell(:);
else
    shocell = shocell(:);
    if length( shocell ) ~= nx
        error('Mistmatch length SHOCELL and SDX HH matrices')
    end
end

% --------------------------------
% Decide whether using DATA or NOT 
% --------------------------------
if nargin < 8 | isempty( in ); 
    in.nrep = 1000; 
end     
if nargin < 9 | isempty( Y )
    flag_data = 0 ;
    nr = 2;
    T = 150 ; 
else
    flag_data=1;
    [T,nvar]=size(Y); 
    nr = 3;
    if nvar ~= nz
        error('Mistmatch matrix of observables and data')
    end
    clear nvar; 
    
    ycov = cov( Y );
    yvar = ( diag( ycov ) )';
    ystd = sqrt( yvar );
    ycorr = corrcoef( Y );
    yacone = accov( Y , 1);

end 

% ====================
%       Top Rows 
% ===================
if nargin < 10 | isempty( info )  
    flag_info = 0; 
else 
    flag_info = 1; 
end 

% ==================================
% Compute moments implied by model %
% ==================================

% --------------------------------
% Check field settings 
% --------------------------------
[ee,in] = ch_field( in , 'nrep' , 1000 ); 
[ee,in] = ch_field( in , 'Ttake'   , T ); 
[ee,in] = ch_field( in , 'Tdrop'   , 200 ); 

Ttake = in.Ttake;

% -----------------------------------------
% CMAT: top cell with Model information 
% ----------------------------------------
if flag_info == 1 
    ncol = max( 6 , nz + 1 );
    
    cmat=emptycell(10,ncol);
    cmat(1,1)={'MOMENTS'};
    cmat(1,4:5)={'SV?',info.sv };
    cmat(2,1:3)={'Authors:','Specification:','Subfolder:'};
    cmat(3,1:3)={info.auth,info.spec,info.subf};
    opath   = extrsubf( info.outpath , cucd );
    % ========================
    % General
    % ========================
    cmat(5,1:2)={opath,'Output PATH'};
    cmat(6,1:2)={info.fmomname,'Output FILE'};
    cmat(7,1:2)={info.tdate,'Date Today'};
    cmat(8,1:2)={info.priornam,'Prior File'};
    % =================================
    % Sample
    % =================================
    cmat(5,4)={'DATA'};
    if isstr( info.demean ) == 1 | iscell( info.demean ) == 1
        cmat(6,4:5)={'Demean?',info.demean }
    else
        cmat(6,4:5)={'Demean?',num2ans( info.demean )};
    end

    cmat(7,4:6)={'Sample',num2str(info.start),num2str(info.finish)};
    cmat(8,4:6)={'Estimation',num2str(info.startest),num2str(info.finishest)};
    cmat(9,1:2)={'Sim Length',num2str( Ttake ) }; 
    cmat(9,4:5)={'Num Rep', num2str( in.nrep ) }; 
    
else 
    
    ncol = max( 4 , nz + 1 );
    cmat=emptycell(3,ncol);
    cmat=emptycell(10,ncol);
    cmat(1,1:2)={'MOMENTS'};
    cmat(3,1:2)={'Sim Length',num2str( Ttake ) }; 
    cmat(3,3:4)={'Num Rep', num2str( in.nrep ) }; 

end 







% --------------------------------
% BASE: Base Cell for output   %
% --------------------------------
addtop = 3; 

base = emptycell( nr + addtop , nz + 1); 
base(3,2:end) = sercell' ; 
base(2,1) = {'============='}; 

rowsim  = addtop+1; 
rowcell = {' SIM MOD';' STA MOD'}; 

if flag_data == 1 
    rowcell =[{' DATA'};rowcell];
    rowsim  = rowsim + 1; 
end 
        
% ============================================
% Compute Simulated moments %
% =============================================

TT=Ttake+in.Tdrop;
nrep = in.nrep; 

covmat = zeros( nz , nz , nrep ); 
acallmat = zeros( nz , nz , nrep ); 
corrmat   = zeros( nz , nz , nrep ); 
% Autocorrelation function (simulated) 
acallmat=zeros(nz,nz,nper+1,nrep);

whandle = waitbar(0,'Computing Simulated Moments'); 

ii=1; 
for ii=1:nrep; 

    waitbar( ii/nrep , whandle ); 
    
    shocks=SDX*randn(nx,TT);

    Ysim=simgensys(shocks,G1,CC,HH,ZZ);
    Ysim=Ysim(:,in.Tdrop+1:end);
    Ysim=Ysim';
    covmat(:,:,ii)=cov(Ysim); 
    corrmat(:,:,ii)=corrcoef(Ysim); 
    aconemat(:,:,ii)=accov(Ysim,1);
    [acallmat(:,:,:,ii)]=accov(Ysim,[0:nper]);
    
end 
close(whandle); 

% ========================================================
% ========================================================
%                     Compile Output                     %
% ========================================================
% ========================================================

% =========================================================
%                   DCELL 2) DATA Moments 
% ==========================================================
if flag_data == 1
    dcell = emptycell( 2, ncol );
    dcell(1,1:2) = {'DATA','MOMENTS'};
    dcell(2,1:2) = {'=========','========'};

    dcell_1 = crtcell( ycorr  , sercell , sercell , {'DATA CORRELATION'} , 1 );
    dcell_2 = crtcell( yacone , sercell , sercell , {'DATA CORR ONE'} , 1 );
    dcell = merge_cells( dcell, dcell_1 , 1 );
    dcell = merge_cells( dcell, dcell_2 , 1 );
    clear dcell_1 dcell_2;

end



% =========================================================
%                   OCELL 2) Simulated Moments 
% ==========================================================
ocell = emptycell( 2, ncol ); 
ocell(1,1:2) = {'SIMULATED','MOMENTS'}; 
ocell(2,1:2) = {'=========','========'}; 

pmat = round( nrep*[0.05;0.95] ); 

vcfull = squeeze( median( covmat , 3 ) );
covmat= sort( covmat , 3 ); 
covmat_1 = squeeze( covmat(:,:,pmat(1) ) );
covmat_2 = squeeze( covmat(:,:,pmat(2) ) );

% ==============================
%      Standard Deviation      % 
% ==============================
sdevf  = sqrt( diag( vcfull ) ); 
sdevf_1 = sqrt( diag( covmat_1) ) ; 
sdevf_2 = sqrt( diag( covmat_2) ) ; 
clear covmat*; 

topcell = emptycell(  6 , nz + 1 ); 
topcell(1,1)={'SIM STD'}; 
topcell(1,2:end) = sercell'; 
topcell(2,2:end) = num2cprec( sdevf' );
topcell(3,1) = {'5%'}; 
topcell(3,2:end) = num2cprec( sdevf_1' );
topcell(4,1) = {'95%'}; 
topcell(4,2:end) = num2cprec( sdevf_2' );

topcell(5,1)={'SIM VARIANCE'}; 
topcell(6,2:end) = num2cprec( (diag( vcfull ))' ); 
clear vcfull; 

if flag_data == 1 
    ocell= merge_cells( dcell , ocell , 1 ); 
    clear dcell; 
    
end 


ocell = merge_cells( ocell, topcell , 1); 
clear topcell; 

stdsim_med = sdevf'; 
clear sdevf*; 

% Autocorrelation 
% ---------------
[aconecell,xtemp] = out_sqmat( aconemat, [] , [] , sercell , [] , 'AUTOCORR ONE' ); 
acsim_med = diag( xtemp.med ); 
ocell = merge_cells( ocell , aconecell , 1);
clear xtemp aconemat aconecell; 

corrcell = out_sqmat( corrmat, [] , [] , sercell , [] , 'CORRELATION' ); 
ocell = merge_cells( ocell , corrcell, 1); 
clear corrmat corrcell; 

% Autocorrelation functions 
% ===========================
subpg=4; 
subc =0; 
pco=0; 
tempx=[0:nper]; 
ii=1; 
for ii=1:nz; 
    temp=reshape(acallmat(ii,ii,:,:),[nper+1 nrep])'; 
    temp=sort(temp); 
    temp_m=mean(temp); 
    temp=temp(pmat,:); 
    empac=squeeze(accov(Y(:,ii),[0:nper]))'; 
    if subc==0 
        figure; 
    end 
    subc=subc+1; 
    
    subplot(2,2,subc); 
    hold;
    plot(tempx,temp_m,'b');
    plot(tempx,temp(1,:),':r',tempx,temp(2,:),':r');
    plot(tempx,empac,'k--');
    title([sercell{ii},' AC mean, [5;95] and data'] );
    hold off;
    axis tight;
    xlim([0 nper]); 
    
    if subc==subpg | ii==nz; 
        subc=0;
        cd(info.outpath);
        saveas(gcf,[sercell{ii},' desc'],'fig');
        print(gcf,'-deps',[sercell{ii},' desc']);
        cd(cucd);
    end 
    
end 

clear sdevf; 

% Model Implied Autocorre
[ee,info]=ch_field(info,'nforc');
[vdcom,vcfull,sdevf,ACzero,ACone,junk,junk,vdech]=vardecom_mode(G1,HH,SDX,ZZ,info.nforc); 

tempind=find( ismember([1:info.nforc]',[1 2 4 8 12]') ); 

% Output for the Variance Decomposition 
vdech=vdech(:,:,tempind); 
temptop=1:info.nforc; 
temptop=temptop(tempind); 

vhcell=emptycell(1,nx+1); 
vhcell(1,1:2)={'VDEC ','Different H'}; 
jj=1; 
for jj=1:size(vdech,3); 
    tempc=crtcell(squeeze(vdech(:,:,jj)),sercell,shocell,['Horizon ',num2str(temptop(jj))] ); 
    vhcell=merge_cells(vhcell,tempc,1); 
end 
clear vdech tempind temptop jj tempc 
    
% =======================================
%        AUTOCORRELATION 
% =======================================
temp_row = strcat( 'AUTOCORR ',rowcell ); 
base( addtop+1:end , 1 ) = temp_row(:); 
clear temp_row; 

base(1,1) = {'AUTOCORRELATIONS ALL'}; 
base( rowsim , 2:end   ) = num2cprec( (acsim_med(:))') ; 
base( rowsim+1 , 2:end ) = num2cprec( diag( ACone )' );  
if flag_data == 1; 
    base( addtop + 1 ,  2:end ) = num2cprec( diag(yacone) )' ; 
end 

ocell=merge_cells(base,ocell,1); 

% =======================================
%        STANDARD DEVIATION 
% =======================================
base(1,1) = {'STD ALL '}; 
temp_row = strcat( 'STD ',rowcell ); 
base( addtop+1:end , 1 ) = temp_row(:); 
clear temp_row; 

base( rowsim , 2:end   ) = num2cprec( stdsim_med' ); 
base( rowsim+1 , 2:end ) = num2cprec( sdevf' ); 
if flag_data == 1; 
    base( addtop + 1 ,  2:end ) = num2cprec( ystd' ); 
end 

ocell = merge_cells( base , ocell , 1 ); 
clear base; 

ocell = merge_cells( cmat , ocell , 1 ); 
clear cmat; 
    
% ==============================
%        MODEL MATRICES 
% ==============================
matcell = emptycell( nz + 2 , nz + 1 ); 
matcell(3:end,1) = sercell; 
matcell(2,2:end) = sercell'; 
tempc = emptycell( nz , nz ); 
indic = find( tril( ones(nz) ) ~= 0 ); 
% AC MATRIX
% ---------
tempc(indic) = num2cprec( ACzero(indic) ); 
matcell(1:2,1) = { 'CORRELATION','MODEL STA' }; 
matcell(3:end,2:end) = tempc ; 
ocell = merge_cells( ocell, matcell , 1);
clear tempc indic ; 

% CORRELATION LAG ONE 
% -------------------
matcell(1:2,1) = { 'AC LAG ONE','MODEL STA' }; 
matcell(3:end,2:end) = num2cprec( ACone ); 
ocell = merge_cells( ocell, matcell , 1); 
clear matcell; 

% VARIANCE DECOMPOSITIONS 
% ========================
vcell = emptycell( nz + 3, nx + 1); 
vcell(1,1:2)={'VARIANCE','DECOMPOSTION'}; 
vcell(2,1) = {'MODEL STA'}; 
vcell(2,2)  ={'Shocks'}; 
vcell(3,2:end) = shocell'; 
vcell(3,1) = {'Series'}; 
vcell(4:end,1) = sercell; 
vcell(4:end,2:end) = num2cprec( vdcom ); 

ocell = merge_cells( ocell, vcell, 1); 
ocell = merge_cells( ocell, vhcell , 1); 

disp('Create cell with model moments'); 
disp('For simulated data and stationary model'); 
